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We present a mechanistic model for a Newtonian fluid called fluid particle dynamics. By analyz- 
ing the concept of "fluid particle" from the point of view of a Voronoi tessellation of a molecular fluid, 
we propose an heuristic derivation of a dissipative particle dynamics algorithm that incorporates 
shear forces between dissipative particles. The inclusion of these non-central shear forces requires 
the consideration of angular velocities of the dissipative particles in order to comply with the con- 
servation of angular momentum. It is shown that the equilibrium statistical mechanics requirement 
that the linear and angular velocity fields are Gaussian is sufficient to construct the random ther- 
mal forces between dissipative particles. The proposed algorithm is very similar in structure to the 
(isothermal) smoothed particle dynamics algorithm. In this way, this work represents a generaliza- 
tion of smoothed particle dynamics that incorporates consistently thermal fluctuations and exact 
angular momentum conservation. It contains also the dissipative particle dynamics algorithm as a 
special case. Finally, the kinetic theory of the dissipative particles is derived and explicit expressions 
of the transport coefficients of the fluid in terms of model parameters are obtained. This allows to 
discuss resolution issues for the model. 



I. INTRODUCTION 



Complex fluid systems such as colloidal or polymeric suspensions, micelles, immiscible mixtures, etc. represent a 
challenge for conventional methods of simulation due to the presence of disparate time scales in their dynamics. There 
is presently a great effort in developing new techniques of simulation that overcome some of the difficulties of micro- 
scopic (molecular dynamics) and macroscopic (numerical solution of continuum equations) conventional techniques. 

On one hand, molecular dynamics captures all the detailed dynamics from times scales of the order of atomic collision 
times to macroscopic hydrodynamic times. However, in order to explore these large macroscopic times the number 
of particles required is enormous. Although large scale molecular dynamics simulations with millions of particles are 
currently performed one realizes that not all the information generated is actually required or even relevant at the 
time scales at which rheological processes in complex fluids take place. 

On the other hand, from a continuum point of view the conventional solution of partial differential equations like the 
Navier-Stokes equation encounters difficulties due to the cumbersome treatment of moving boundary conditions to be 
imposed in a system as, for example, a colloidal suspension. These problems can be alleviated by the use of Lagrangian 
descriptions in which the discretizing grid moves according to the flow. A particularly exciting development has been 
the technique of smoothed particle dynamics (SPD) which is essentially a discretization by weight functions that 
transforms the partial differential equations of continuum mechanics into ordinary differential equations [jl],^) . These 
equations can be further interpreted as the equations of motion for a set of particles interacting with prescribed laws 
of force. The technique thus allow to solve PDE's with molecular dynamics codes. Another advantage of a Lagrangian 
description relies on the fact that no expensive recalculations of the mesh are required as the dynamics takes care 
of it. Smoothed particle dynamics has been used for simulating astrophysical flows with non- viscous terms (this was 
the original aim of the technique at the early 70's 0), and very recently in the study of viscous H|| and thermal 
flows |^J^] in simple geometries. It has not been applied to the study of complex fluids and this may be due, in 
part, to the fact that there is no easy implementation of Lagrangian fluctuating hydrodynamics |Q with SPD. Such 
implementation would be highly desirable in order to study the Brownian realm in which many of the processes in 
complex fluids take place. It must be noted that the particulate nature of the algorithm in SPD produces fluctuations 
which, from a computational point of view, are regarded as numerical noise. It is not clear that in the presence of 
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viscous terms this noise satisfies the appropriate fluctuation-dissipation theorem. A second problem with SPD is that 
for viscous problems the non-central nature of the viscous shear force breaks the conservation of angular momentum, 
even though the initial continuum equations are perfectly isotropic and conserve angular momentum. 

In between microscopic and macroscopic descriptions, mesoscopic levels of description are gaining attention in order 
to address flow problems in complex fluids and/or geometries. Lattice gas automata [^[^), lattice Boltzmann automata 
fiof or the direct simulation Monte Carlo method for dilute gases 1 11 have been useful tools in studying hydrodynamic 
problems in complex geometries. For the case of colloidal suspensions, lattice Boltzmann techniques represent a serious 
competitor to Brownian dynamics fl^l or Stokesian dynamics |l3f| in that the computational cost scales linearly with 
the number of colloidal particles whereas, as a consequence of the l ong ranged hydrodynamic interactions, it increases 
with the cube of the number of particles in the latter techniques |10| . A drawback of the lattice approaches is that 
the dynamics is constrained by the lattice. This makes the consideration of boundary conditions on shaped bodies 
cumbersome. 

In the same spirit of looking at mesoscopic descriptions, a very appealing idea was introduced by Hoogerbrugge 
and Koelman [fujjl^] in which a coarse grained description of the solvent fluid in terms of dissipative particles was 
devised. The technique was coined dissipative particle dynamics (DPD) and it is an off-lattice technique that does not 
suffer from the above mentioned drawbacks of lattice gas and lattice Boltzmann simulations. DPD consists essentially 
on a molecular dynamics simulation in which the force between particles has, in addition to a conservative part, 
a dissipative part represented as a Brownian dashpot. This Brownian dashpot damps out the relative approaching 
velocity between particles and introduces a noise term that keeps the system in thermal agitation. The dissipative 
particles are understood as "droplets" or cluster of molecules that interact with each other conserving the total 
momentum of the system |l4|,[l6|]. This global conservation law has its local counterpart in the form of a balance 
equation for the momentum density, and the dissipative particles behave hydrodynamically in the low wave number 
and frequency regime. 

DPD has received a substantial theoretical support. It has been shown that the original DPD algorithm of Ref. |l4|] 
has associated, under a slight modification, a Fokker-Planck equation with Gibbs equilibrium states [17]]. The extension 
to multicomponent systems has also been considered |l8f . A first principles derivation of DPD for a harmonic chain 
has been presented in ]l9[ |. The macroscopic hydrodynamic equations have been obtained with projection operator 
techniques J2^|. A very important further step has been the formulation of the kinetic theory for DPD by Marsh, 
Backx, and Ernst pl[ , which allows to relate the transport coefficients in the hydrodynamic equations with the DPD 
model parameters. Finally, the effect of finite time steps on the equilibrium state of the system has been considered in 
[ 22 1 . DPD has been since applied to the study of colloidal suspensions |l5|]2"3| , |2"4| , porous flow [jlij, polymer suspensions 
[25 j, and multicomponent flows pq| . 

In this work we provide a more precise meaning to the concept of "droplet" or "fluid particle" from a Voronoi 
tessellation of physical space. This conceptual framework allows to model the different processes that intervene in 
the interaction between fluid particles or mesoscopic clusters of atoms of the fluid. The outcome is a generalization of 
the algorithm of DPD that includes shearing forces between the fluid particles. These forces are non-central and do 
not conserve total angular momentum. This enforces the inclusion in the model of a spin variable with a well-defined 
physical meaning for the fluid particles. In this way angular momentum conservation is restored. We also investigate 
the structure of the random forces that must be included in order to recover a Gaussian distribution of linear and 
angular velocities for the fluid particles (note that the equilibrium fluctuations of the hydrodynamic velocity field are 
Gaussian). The structure of the random forces is postulated after analogy with the structure of the random stress 
tensor in terms of the Wiener process in the fluctuating hydrodynamics theory |^8| . 

We show that the form of the equations of this fluid particle model at zero temperature (when fluctuations are 
absent) and with no angular variables is identical to the form of the equations obtained in a simple version of smoothed 
particle dynamics as applied to fluid systems. In this sense, this work can be regarded as a generalization of SPD 
that includes fluctuations consistent with the principles of statistical mechanics and that conserves exactly the total 
angular momentum of the system. In other words, the obtained fluid particle algorithm may be viewed as a Lagrangian 
discretization of the equations of (isothermal) fluctuating hydrodynamics. 

The paper is structured as follows. Section II considers the definition of a fluid particle in terms of the Voronoi 
tessellation and this serves to motivate the type of forces and torques between fluid particles introduced in Section 
III. Section IV presents the equivalent Fokker-Planck equation to the equations of motion. This allows to establish 
requirements on the model parameters in order to have a proper equilibrium distribution. A summary of the model 
is presented in Section V. Section VI contains the kinetic theory of the model in the simple case when conservative 
forces are absent. The transport coefficients are computed and this permits to discuss its dependence on the number 
density of fluid particles in Section VII. A final discussion and conclusions is presented in the last section. 
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II. FLUID PARTICLES THROUGH VORONOI TESSELLATION 



In an attempt to better understand the physical meaning of DPD, we have devised a coarse graining procedure 



for a molecular dynamics simulation of point particles (atoms) interacting through continuous potentials |16|. The 
coarsening is performed through the Voronoi tessellation that allows to divide physical space into a set of non- 
overlapping cells that cover all the space in a well-defined manner. Given a discrete set of points (that can be 
distributed at random) the Voronoi tessellation assigns to each point (called "cell center"), that region of space that 
surrounds it and that is nearer to this point than to any other point of the set. With this tessellation the atoms of 
the molecular dynamics simulation are distributed into clusters around the centers of the Voronoi cells. The practical 
way to perform the Voronoi tessellation in the simulation is by computing the distance of a given atom to all the 
center cells and assign that atom to the nearest center. Subsequently, the Voronoi cells are set in motion according 
to the velocity and acceleration of the center of mass of the particles that are within the cell. In this way, the cells 
capture the motion of the fluid at mesoscopic scales. 

The Voronoi cells are a well defined realization of what is loosely regarded in fluid mechanics textbooks as "fluid 
particles". We would like to know how these fluid particles move, that is, which explicit law of force between fluid 
particles would reproduce the actual motion of the Voronoi cells observed in the simulations. It is apparent that the 
number of cells is much smaller than the number of atoms in the molecular dynamics simulation, and therefore if 
one knows how the clusters move, one can try to simulate the clusters and capture the mesoscopic behavior of the 
underlying liquid with much less computational effort. 

For sufficiently large clusters, that is, when the typical distance A mes between cell centers is much larger than the 
typical distance \ m ic between atoms we expect that the clusters move hydrodynamically. More precisely, we consider 
the conserved densities 

M r 

& ~ V T 

where M r ,P r ,£ r are the instantaneous total mass, momentum and energy, respectively, of the system of particles 
that happen to be within the Voronoi cell centered at r and V r is the volume of the cell. One also introduces 

v r = ^ (2) 

Pr 

which is the instantaneous velocity of the center of mass of the system of particles within the Voronoi cell at r. 

If, 1) the cells are large enough for the system of particles that are within it to be considered as a thermodynamic 
system, and 2) the variations of the conserved quantities from neighbor cells are small, then the variables (|l|) obey the 
equations of fluctuating hydrodynamics @] (see (28) for the non- linear case). The conserved quantities (0) are subject 
to fluctuations because the atoms can enter and go out from the Voronoi cells due to their thermal agitation. The size 
of fluctuations, that is, the noise amplitude appearing in the equations of fluctuating hydrodynamics is proportional 
to the square root of the inverse of the volume of the cell, in accordance with the 1 / ^/N dependence of fluctuations 
in equilibrium ensemble theory pgf| . Therefore, depending on the "resolution" (the number of Voronoi cells per unit 
volume) used to describe the system, the amplitude of the noise term in the hydrodynamic equations that govern the 
instantaneous values of the conserved variables will be different. 

Now, one is faced with two possible routes in order to simulate the dynamics of clusters. The first route is to 
consider the conserved discrete variables (|l]) as the state variables and update them according to some discretized 
version of the equations of hydrodynamics. This poses some subtle problems regarding the formulation of fluctuating 
hydrodynamics in a moving mesh, in particular with the treatment of the l/\/Vr singularity. The second route, which 
is the one we follow in this paper, is to postulate the laws of force between cells. Despite the strong assumptions 
made to model these forces, the final expressions satisfy symmetry requirements that ensure that the behavior of the 
clusters will be, on average, that of real clusters. 

A final word on the angular momentum is in order. We have not included in the above set of conserved variables 
(|l|) the angular momentum density defined by 

g 

J r = r x g r + — (3) 
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where we have decomposed the angular momentum of the system of particles in cell r as the sum of the angular 
momentum of the center of mass with respect to the origin plus the intrinsic angular momentum S r of the particles 
of the cell with respect to the center of mass of the cell. We could, in principle, compute the inertia tensor I r for the 
system of particles in cell r and define the angular velocity oj r through 

= I- X S r (4) 

The reason why the angular momentum density is usually not considered in the derivation of Newtonian hydrody- 
namics is because the second contribution in (||) vanishes in the "continuum limit" . This can be seen by noting that 
S r must scale as a typical size of the cell. In the continuum limit this typical distance goes to zero and there is no 
intrinsic angular momentum contribution. The situation here is different from the case of molecular fluids with spin 
p7f , where the rotation of the molecules themselves originates an angular momentum that does not scale with the 
size of the cells and produces a finite value in the continuum limit. 



III. MODELING THE FORCES AND TORQUES BETWEEN FLUID PARTICLES 

In this section we formulate the fluid particle model under a set of simplifying hypothesis. In the real clusters, the 
mass is a fluctuating quantity as particles can enter and go out from the Voronoi cell. Also, the shape of the cells 
changes as the cells move. However, the first assumption is that all clusters are identical, having a fixed mass m and 
fixed isotropic inertia tensor of moment of inertia /. We assume that the state of the cluster system is completely 
characterized by the positions r^, the velocities of the center of mass and the angular velocities uji. Note that we 
do not include any internal energy variable and therefore, the resulting algorithm will not capture appropriately the 
thermal effects that occur in real fluids. This may be a minor problem when one is interested only in rheological 
properties. A generalization of the model including energy conservation has been recently proposed independently in 
Refs. [0 (29). 

The next step is to specify the forces and torques that are responsible for changing the values of the linear and 
angular velocities of the clusters. We model the forces between two clusters by considering several heuristic arguments 
about how one expects that the actual Voronoi clusters interact with each other. In this respect we make first a strong 
pair-wise additivity assumption. In the real system one expects that the force between two clusters (that is between 
all the atoms of the first cluster that are interacting with the atoms of the second cluster) will depend in general not 
only on the state variables of these two clusters but also on the configuration of other neighboring clusters. For the 
sake of simplicity, though, we neglect this collective effect and assume that the force between two clusters depends 
only on the position and velocities of these two clusters. 

The equations of motion are therefore 

Vi = V s ; 

V 4 = — Yl F ij 

m ^ 

Wi = j X] N « ( 5 ) 

where Fy ,Ny are the force and torque that cluster j exerts on cluster i. We require that the forces satisfy Newton's 
third law, Fy = — F^j, in such a way that the total linear momentum P = mv; is a dynamical invariant, P = 0. 
In addition, we assume that the torques in (|^) are given by 

1 
2 

and one checks immediately that the total angular momentum 



N « = T7 r <. x F . ( 6 ) 



J = X Tj X Pi + Iu>i (7) 

i 

is conserved exactly, J = 0. 

We will model the force between clusters according to 

F id = Ffj + F?j + F* + F tJ (8) 

The first three contributions are deterministic forces whereas the last one is random. Let us discuss them separately. 
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A. Deterministic forces 



The first contribution Ffj to the force is assumed to arise from a conservative potential V(r) that depends on the 
separation distance between clusters. In Ref. |l6| we have argued that a plausible definition of this potential is through 
the logarithm of the radial distribution function of cluster centers. The resulting soft potential has a bell-shaped form 
and has the virtue that when used as the potential between clusters in a MD simulation, it reproduces consistently 
the radial distribution function of the real clusters (as it has been checked through an actual MD simulation). It 
therefore captures the static or equilibrium properties of the system of clusters. The physical interpretation of this 
force is that it provides the excluded volume effect of each cluster. The center of the cluster (and its center of mass) 
is usually located "in the middle" of the cell and therefore it is not very probable that two cluster centers are closer 
to each other than the typical size of the cell. 

It is clear, though, that this conservative contribution cannot be the only contribution to the force because it does 
not capture friction effects between clusters. These friction effects will depend on the velocities between clusters and 
will give rise to dissipative processes. The second contribution in (j^) is a friction force that depends on the relative 
translational velocities of the clusters i, j with positions r^r^ and velocities Vj,Vj in the following way 

F£ = - 7 mM T (r y -)-v. y - (9) 

where Vjj = Vj — Vj is the relative velocity and the dimensionless matrix 'M T (r i j) is the most general matrix that can 
be constructed out of the vector = r , — , this is 

M T (ry ) ee A(ry)l + Bfojjeyey (10) 

where 1 is the unit matrix, e,j = ry/r^j is the unit vector joining the particles, = |ry| and the dimensionless 
functions A(r) and B{r) provide the range of the force. The friction coefficient 7 has been introduced as an overall 
factor for convenience and has dimensions of inverse of time. The first contribution to the dissipative force (j^) is in 
the direction of the relative velocities and tends to damp out the difference between the velocities. It is a shearing 
force which is non-central. The second contribution is directed along the joining line of the particles and damps out 
the relative approaching motion of the particles. The dissipative force in the original algorithm of Hoogerbrugge and 
Koelman is obtained with A(r) = |Q. Note that the form of the force Ffj is the more general expression for a 
vector that depends on and is linear in the relative velocities. 

We now discuss the effects of rotation in the dissipative force. Let us assume for a moment that the clusters i and j 
were spheres of radius /2 in contact and spinning with angular velocities Wj, cjj with no translational velocities. We 
would have a relative velocity at the "surface" of the spheres equal to x (u>i and it is plausible to associate 

a friction force between the spheres proportional (in matrix sense) to this relative velocity. Therefore, the rotational 
contribution to the dissipative force is of the form 

F* = - 7 mM s ( rij ). x [ui + uj]) (11) 

Again, the dimensionless matrix depends only on the vector and therefore it must have the form 

M R (r tJ ) = C{r l3 )l + DfoOeyea (12) 

where C(r), D(r) are scalar functions. The first part of the matrix gives rise to a friction force proportional to the 
relative velocity at the "surface" of the spheres. The effect of this force is twofold. On one hand, the spinning of a 
particle causes translational motion onto the neigbouring particles. On the other, it also causes rotational motion in 
such a way that two neigbouring particles prefer to be with opposite angular velocities (in a sort of "engaging" effect) . 
The presence of a third particle frustrates the spinning of both particles and, therefore, the global effect of this force 
is to damp out to zero the angular velocities of the particles. The second contribution to the force (|l^) is actually 
zero because the cross product is perpendicular to e^ . We retain this term just to maintain the analogy between both 
matrices in (^|) and M T in (|TT|). Finally, we use the same value for 7 in (||) and in (jlj) because any difference 
can be taken into account through the functions A(r), B(r), C(r), D(r). 

If we use, instead of the polar vector representation for the angular velocity, the antisymmetric tensor representation 
tu xy = —TjJy X = lo z (cyclic), we can write the force in the form 

Fg - - im m R {v l3 ) ■ { Wi + w s y r -f (13) 

which shows explicitly the vectorial nature of the force (this is, F^ transforms under rotations as a vector). 



5 



B. Random forces 



The first three contributions Ff^ , Ffj , to the force between clusters in (||) are deterministic while the last one is 
stochastic. The reason why we introduce a random force is because it is well-known that whenever a coarse-graining 
procedure is performed, dissipation and noise arise and both are related through a fluctuation-dissipation theorem. 
After discussing the form of the dissipative forces we will now consider the form that the random force should have. 

Inspired by the tensorial structure of the random forces that appear in the fluctuating hydrodynamics theory 
p8| , we expect that the dissipation due to shear has associated a traceless symmetric random matrix and that the 
dissipation due to compressions has associated a diagonal trace matrix. By symmetry reasons, we expect that the 
noise associated to rotational dissipation will involve an antisymmetric matrix. Therefore, we postulate the following 
velocity independent random force 

Fydt = am (l(r y )dWj + /,': r, ; r ,/\V, ; 1 + C^)^ ) e, ; ( 1.4) 



where A(r), B(r), C(r) are scalar functions, a is a parameter governing the overall noise amplitude, and we introduce 
the following symmetric, antisymmetric and traceless symmetric random matrices 



dW^ = I[dW^-dW^] 



rfW - = dW?- - -tr[dW?-]l (15) 



D 

The overline in a matrix denotes its traceless part. Here, D is the physical dimension of space and dW^ is a matrix 
of independent Wiener increments which is assumed to be symmetric under particle interchange 

dWg = dWg (16) 

This symmetry will ensure momentum conservation because F^ = — F^. The matrix <iW^ is an infinitesimal of 
order 1/2, and this is summarized in the Ito mnemotcchnical rule 

From this stochastic property, one derives straightforwardly the following rules from the different parts in ( |14| ) 
tr[dW«']tr[dW#/] = [8^S Vf + 6^6^] Ddt 

1 1 

dt 



dW%? dW%, = [SySyj, + 8 ir 8 jV ] S^S^dt (17) 



dW^r dWj-" = [SijSvj, +6 ij >8 ji > 



1 



2 

tridWi^dWjj, = 
tvidWu^dWfj, = 

dwlr'dwf^' = o (is) 

These expressions show that the traceless symmetric, the trace and the antisymmetric parts are independent stochastic 
processes. The apparently complex structure of the random force ( [lij ) is required in order to be consistent with the 
tensor structure of the dissipative friction forces (^J) and (|Tl|). This will become apparent when considering the 
associated Fokker-Planck equation in the next section and requiring that it has a proper equilibrium ensemble. 

Despite of the heuristic arguments and strong assumptions made in order to model the force Fy between clusters, 
we note that this force is the most general force that can be constructed out of the vectors r,, r,-, Vj, Vj, ojj,, ojj and 
that satisfies the following properties: 

1. It is invariant under translational and Galilean transformations and transforms as a vector under rotations. 

2. It is linear in the linear and angular velocities. This linearity is required in order to be consistent with the 
Gaussian distribution of velocities at equilibrium, as we will show later. 

3. It satisfies Newton's third law F^- = — Fji and, therefore, the total linear momentum will be a conserved quantity 
of the system. 
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IV. FOKKER-PLANCK EQUATION AND EQUILIBRIUM STATE 



The equations of motion are Langevin equations which in the form of proper stochastic differential equations 
(SDE) become 



dxi = Widt 

m — ' L 



i' i' 



(19) 



where we have introduced 



</v„- = ifii/t// = a ( _!(/•;,- ></W; v + ij( / „v)-^tr[,/W,/l + C(r,,ulW;), j-e 



(20) 



In p rinciple, one should specify which stochastic interpretation (Ito or Stratonovich) must be used in these equations 
Nevertheless, both interpretations produce the same answers because the random forces are velocity independent. 
Associated to the SDE (|19|) there exists a mathematically equivalent Fokker-Planck equation (FPE). The FPE 
governs the distribution function p(r, v, to; t) that gives the probability density that the N clusters of the system have 
specified values for the positions, velocities and angular velocities. We show in the Appendix that the FPE is given 
by 



d t p(r, v, to; t) = [L c + L T + L R ] p{r, v, u; t) 



(21) 



The operator L is the usual Liouville operator of a Hamiltonian system interacting with conservative forces F , this 



is, 



E^ + E^t 



(22) 



The operators L T , L R are given by 



with 



T T _ -pT 



2 J 



(23) 



(9 



d 



ma" 
r . . 

I 2 y 



Here, the matrix Ty is given by 



A 2 { nj ) + C 2 {r i3 ) 









1 + 


[G- 





d d 
duji dujj 



i 2 (^-) + ^ 2 (ry)-^ 2 (ry) 



(24) 



(25) 



The steady state solution of equation (^lj), dtp — 0, gives the equilibrium distribution p eq . We now consider the 
conditions under which the steady state solution is the Gibbs canonical ensemble: 



p eq (r,v,ij) = |exp{-F(r,«,cj)/fc s T} 



1 , / \ - TOn I 



+ K(r) /fc B T} 



(26) 
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where H is the Hamiltonian of the system, V is the potential function that gives rise to the conservative forces F c , 
ks is Boltzmann's constant, T is the equilibrium temperature and Z is the normalizing partition function. We note 
that the velocity and angular velocity fields are Gaussian variables at equilibrium and, therefore, one expects that the 
distribution function of the discrete values of theses fields is also Gaussian. 

The canonical ensemble is the equilibrium solution for the conservative system, i.e. L c p eq = 0. If in addition the 
following equations are satisfied 

^Hjp eq = 

L*p eq = (27) 

then we will have Lp eq — and the Gibbs equilibrium ensemble will be the unique stationary solution of the dynamics. 
Eqns. (p7|) will be satisfied if 



a 2 m 



(28) 



' 2k B T 

which is a detailed balance condition, and also 

M>y) = M T (r y ) = Ty (29) 
This is the fluctuation-dissipation theorem for the fluid particle model. We observe, therefore, that the initial hy- 



pothesis for the tensorial structure of the dissipative (|9|) , ( |13| ) and random (|1J) forces was correct and consistent with 
equilibrium statistical mechanics. 

A final word about an H-theorem is in order. It has been shown in Ref. that the original DPD algorithm has an 
H-thcorem that ensures that the equilibrium ensemble is the final solution for whatever initial condition selected. In 
the model presented in this paper there is also a functional of p(z) that is a Lyapunov functional. It is not necessary to 
prove again that an H-theorem exists for the fluid particle model, because a general H-theorem exists for any Fokker- 
Planck equation |32] ]. The only condition is that the diffusion matrix accompanying the second derivative terms of 
the FPE is positive (semi) definite. However, in the model presented in this paper the diffusion matrix is positive 
semidefinite by construction, because the FPE has been obtained from a SDE. The diffusion matrix is obtained from 
the product of two identical matrices. Therefore, its eigenvalues are the square of the eigenvalues of these matrices 
and are necessarily positive (or zero). 



V. SUMMARY OF THE FLUID PARTICLE MODEL 

In this section and for the sake of clarity we collect the results presented so far. The fluid particle model is defined 
by N identical particles of mass m and moment of inertia I. The state of the system is characterized by the positions 
velocities Vj, and angular velocities cjj of each particle. The forces and torques on the particles are given by 



N, = -EfxF tJ (30) 

3 

where the force that particle j exerts on particle i is given by 

F^Fg + F^. + F^ + F^ (31) 
The conservative (C), transnational (T), rotational (i?) and random (~) contributions are given by 



I'l; = -V'{n )e K 
FS = -7mT«.v 



Fg = - 7 mT ir (| x (Wi+wj) 



Fijdt = {2k B T im ) 1/2 ( A( ri] )dW- 3 + B(r y )-^tr[dWy]l + rir^u/W,] ) -ey (32) 
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The random bits are defined in ( |l5| ) and its stochastic properties are given in (|T^). Here, the matrix T is given by 



A( rij )l 



B(jij )ejj-ey 



(33) 



where 



A{r) = -[A\r)+C 2 (r) 
B{r) = \[A*{r)-C\r) 



1 

D 



B 2 {r)-A 2 (r) 



(34) 



The model is thus specified by providing the scalar functions V(r), A(r), B(r), C(r). We note that the case A(r) = 
C(r) = corresponds to the original DPD algorithm of Hoogerbrugge and Koelman jl4|,|l7]]. In this case, the random 
force is given in terms of a single random number (the trace), the forces are central and the torques vanish, rendering 
the spin variables unnecessary. Note that there is some freedom in selecting the functions A(r) , B (r) , C (r) and it 
might be convenient to take A(r) or C(r) equal to zero in order to compute only four of seven random numbers in 
each step of a simulation. 

The model presented in the language of SDE is more appropriate for its direct use in simulations. For theoretical 
analysis it is more convenient to use the corresponding FPE which is given by 



d t p(r, v, to; t) = [L c + L T + L R ] p(r, v, u; t) 



(35) 



where 



Ff _d_ 

m <9vi 



£0 ~2>5J 



L = — 



(36) 



Here, the vector operators are given by 



7T, 



k B T 


' d 


d ' 




m 


dvi 







(|x h+Wj ]) + f f 



d 



d 



(37) 



were the last equality defines the two vector operators V-y, W,. 



VI. KINETIC THEORY 



One would like to predict the macroscopic behavior of the fluid particle model and, in particular, check that this 
behavior conforms to the laws of hydrodynamics (as expected from symmetry considerations) and predict the value 
of the transport coefficients in terms of model parameters. The global conservation laws of mass, linear and angular 
momentum in the fluid particle model have a local counterpart in the form of balance equations. Our aim is to 
formulate these equations of transport within a kinetic theory approach, as has been done by Marsh, Backx, and 
Ernst recently for the case of the original DPD model in Ref. [pi]. A derivation of the hydrodynamic equations 
with a projection operator technique for the original DPD algorithm was presented in Ref. |2Cf| . The projector used 
was the Mori projector |}3| and the resulting equations were the linearized equations of hydrodynamics. By using a 
time- dependent projector one can obtain the non-linear equations of hydrodynamics with the transport coefficients 
expressed in terms of Green-Kubo formulae [p4| . Although explicit calculations can be performed of these Green- 
Kubo formulae under certain approximations |35j , we adopt in this paper the approach of kinetic theory, allowing for 
a straightforward comparison with the results of Ref. [plj . 
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A. General rate of change equation 



The starting point is the formulation of the general rate of change equation for an arbitrary function G(z) where z 
is a shorthand for the set of all positions, velocities and angular velocities of the N particles of the fluid. By using 
the Fokker-Planck equation ( |35| ) , we can write 

d t (G) = J dzG{z)d tP (z-t) 

= J dzG{z) [L c +L T + L R ] p(z; t) 

= J dzp(z; t) [L c+ + L T+ + L B+ ] G(z) (38) 
where an integration by parts is performed and the adjoint operators are defined by 



L 



c+ 



d Ff d 



m <9v 
T 



d 



LS+ = _™ 7E( s xT , i . ( v lj+ w s) ).£- 



Here, the vector operators are given by 



W2 



d 



, k B T 

H 

m 

(Y x [w» + 



<9v 



_9_ 



(39) 



(40) 



to be compared with (| 



B. Balance equations 

The conserved density fields are expected to behave hydrodynamically. The conserved density fields are the mass 
density p(r, t) — mn(r,t), where n(r, t) is the number density field; the momentum density p(r, i)u(r, £), where 
u(r,t) is the velocity field; and the total angular momentum density field J(r, t) = L(r, t) + S(r, t) where L(r, t) — 
r x p(r, t)u(r, t) is the macroscopic angular momentum density and S(r, t) — In(r,t)Cl(r,t) is the intrinsic angular 
momentum density or spin density. Here £l(r, t) is the angular velocity field. The number density, the velocity and 
angular velocity fields are defined by 

n(r,i) = <£*(r -ri)) 

i 

n(r,t)u(r,i) = v;<5(r - r;)) 

i 

n(r, t)fi(r, t) = (£ ^6(r - r 4 )> (41) 

i 

By applying (pq) to the mass and momentum densities (|4l]) we obtain the set of balance equations 

dtp = -Vpu 
d t pvL = — V [puu + n] 

(42) 

where the total stress tensor II = II A + Yl c + U D and the kinetic, conservative and dissipative contributions to the 
stress tensor are 
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n 



ir 



K 



(m^(v 4 - u(r,t))(vj - u(r. t))S(r - r<)) 

i 



n 



D 



dA<5(r -rj - Ary)) 



(43) 



Here, = vy + x [cjj + ujj]/2 is the relative velocity at the "surface of contact" of two identical spheres separated 
a distance TVj . In deriving these equations we have used the identity 



f 1 d 

S(r - r<) - <5(r - r 3 -) = y dA— <5(r - r> + Ar.^ ) 

= -V-r y / d\S(r - r, : + Ar y ) 
Jo 



(44) 



We note that the kinetic and conservative parts to the stress tensor are symmetric tensors but the dissipative part is 
not and therefore we must be careful with the ordering of the indices. In Cartesian components we understand the 
momentum balance equation (f42|) as follows (summation over repeated indices is implied) 



and the dissipative stress tensor has the form 



n 



D 



- 7 m< 



\ E < T -;x / dXS ( r Ar «)) 



Concerning the angular velocity field, by using again Eqn. fl38| ) on the definition (|4lj) we obtain 
d t nn = -V<£ W(r - r,)) + y 7 <£ X ''''-<•'.) 6(r ~ r ^ 



(45) 



(46) 



(47) 



Note that the rate of change of the spin S = Jnfl cannot be expressed entirely as the gradient of a flux. This is 
a reflection that the intrinsic angular momentum S is not a conserved quantity. In the same way, the macroscopic 
angular momentum L = r x pu is not conserved either, as can be appreciated by taking the cross product of the 
momentum balance equation ( |42|) with the position vector r, this is 



<9jL = -r x V(puu + II) 

= -V(Lu + r x n) + 2n A 



(48) 



where H A is the antisymmetric part of the stress tensor (expressed here as an axial vector, this is U Aa — ^e a,lv Tl fJ ' 1 ' 
where e a ^ v is the Levi-Civita symbol). If the stress tensor is symmetric (i.e. its antisymmetric part is zero), the 
macroscopic angular momentum is conserved. In the fluid particle model the non-central nature of the forces implies 
that the antisymmetric part of the stress tensor is not zero. Actually, it is given by (as an axial vector) 



2U A = - 7 m< ]T {^-f x T,,.g,/) f d\S(r - r* + Ar y )) 



(49) 



If we add the last term of (M7j) with the last term of (f48f) which is ( 49 ) , we obtain 



2U A + 7 m<£ x T irSij ) 5(r r. t )> = 7 m<£ x T ir g, 



S(r - Vi) 



I d\S(r — Ti + Xri 
Jo 



7 mV( (y x T u -8ij) r « dX d\'6(r - r, + XX' m)) (50) 



where we have used again (0). Therefore, the total angular momentum density J = L + S satisfies a balance equation 
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d t 3 = -V [Jv + r x n + *] (51) 

where 



(52) 



C. Balance equations in terms of distribution functions 

It is convenient to express the quantities appearing in the balance equations in terms of the single particle and pair 
distribution functions, defined as 

/(x,t) = /(r,v, W ,i) = <I>(x-Xi)> 

i 

/< 2 )(x,x',i) = (]T ^x-x^x'-x,)) (53) 



The number density, the velocity and angular velocities in ( |41| ) are the first moments of the single particle distribution 
function, 



n(r,t) — j dvdujf(r,-v 7 u: : t) 
n(r, <)u(r, t) — J dvdLuvf(r,v,oj,t) 



n(r,f)n(r,t) = j dwdoju)f{v,v,uj,t) (54) 
Next, by using that for an arbitrary function G 

Cy^j G(rij,Vi,vj,uJi,ujj)6(r — Ti + \rij)) = J d-vdV diadio'dllGCR, v, v', u, uj') 

x /( 2 ) (r + AR, v, u, r + (A - 1)R, v', uj') (55) 

which, for A = becomes 

(^~^ Gjrjj, Vj, Vj,u}j,ujj)5(r — r^)) = J d~vdv'duiduj'dRG(R,, v,v',cj,a/) 

x /( 2 )(r,v,w,r-R,v',w') (56) 
we can write the different contributions (|4^) to the stress tensor in terms of the distribution functions 

n A = J dvdujmiy — u)(v — u)/(r, v, u, t) 

n c = y dvdujdv'duj' J dR^F c (R)7 (2) (r,v,w,r',v',w') 

U D = -jm J dvdwdv'du' f dRyT(R)-g/ (2) (r, v, w,r', v',a/) (57) 

where g = [v — v' + ^ x [w + u;']] . We have introduced in these expressions the spatially averaged pair distribution 
function 

7 (2) (r, v, uj, r', v', of) = [ d\f™ (r + AR, v, w, r + (A - 1)R, v', w') (58) 
Jo 

In terms of the distribution functions the terms of the right hand side of ([l7|) can be written as 
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Finally, 



where 



v l w J ;(5(r - r,;)) = J dvdu>vu)f(r,v,u,t) 

= J dvduj (v — u)(oj — fl)f(r, v, ui, t) — nflu 
+ U J dvdujLjf(r,\-,Lu,t) + fl J dvdojvf(r,v,uj,t) 
C£, (Y x Tygy) S(r-n)) = J dvdujdVdu' J dR f~ x T(R)-g) 

x /( 2 )(r,v,c,r-R,v',c') 

* = 7 m y dvdudv'duj' J dR ^ x T(R)-g^ R7 (2) (r,v,w,r',v',w') 

7 (2) (r,v,w,r',v',w') = / dA / dA'/ (2) (r + AA'R, v, w, r + (AA' - 1)R, v', w') 
Jo Jo 



(59) 



(60) 



(61) 



D. Fokker-Planck-Boltzmann equation 

The Fokker-Planck-Boltzmann equation (FPBE) is an approximate kinetic equation for the single particle distri- 
bution function /(x, t). The FPBE is obtained by applying the general rate of change equation d3§| ) to /(x, t). After 
some algebra one arrives at 

dtf + v- V/ = /" dRdv'dw'd- [F C (R) + 7 T(R) -g] /< 2 )(r, v, w, r - R, v', «') 

+ —7 / dRdvWd-T(R)-d/ (2) (r,v,w,r-R,v>') (62) 

TO J 

where we have defined the operator 



_ d to /R d 



In obtaining (|62[), we have inserted at some point the identity 

1 = y dRdVduj'5{r R r,-)£(v' - Vj)<J(<y - w,-) 



(63) 



(64) 



Equation ( |62j ) is not a closed equation for f(r,v,u>,t) because the pair function /( 2 )(x,x',i) appears. Nevertheless 
it can be closed approximately by using the molecular chaos assumption. In what follows we will assume that the 
friction 7 is so large to allow for a neglection of the conservative forces F c H]. This simplifies considerably the 
calculations in the next section. The molecular chaos assumption in the absence of conservative forces becomes 



/< 2 >(x,x',*)«/(x,t)/(x',i) 



(65) 



The final closed Fokker-Planck-Boltzmann equation for the distribution function is, after using the molecular chaos 
assumption (p5[) 



dtf + v-Vf = /[/]= 7 y dRdvW/(r-R,v',^')3-T(R). 
which is an integro-differential non-linear equation. 



gH <? 



/(r,v,cj) 



(66) 
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E. Chapman-Enskog solution of the FPBE 



Our aim is to solve the nonlinear FPBE (66) by using the perturbative method of Chapman and Enskog. The 
method is valid for situations in which the macroscopic conserved fields are slowly varying in typical molecular length 
scales. In these situations, the distribution function decays in a very short kinetic time (short compared to typical times 
of evolution of the conserved field) towards the so called normal solution where the distribution function /(v, oj\a(r,t)) 
depends on space and time only through the first few moments a(r, t) During this last hydrodynamic stage, the 
solution can be obtained perturbatively as an expansion in gradients, this is /(v, w|u, fi) = fo + f\ + 0(V 2 ) where fo 
is of zeroth order in gradients and f\ is of first order in gradients. By substitution of this expansion into the FPBE 
(p6|) one obtains 



9t/o + 9t/i + v-V/ = /[/o] + ( ^) /i + G(V 2 ) «i7i 

fo 



df 



By analogy with the conventional kinetic theory and also with the kinetic theory for DPD in Ref . , we expect that 
the lowest order contribution fo is given by the local equilibrium form for the distribution function. In the presence 
of spin variables it takes the form 

fo(r,v,u,t)=n(r,t)( ™ \ ' exp { —( w - u f\( — _ — ^ ' exp ( —(uj-flfX (68) 

J0K ' ' ' ; V ' ' \2Trk B T J P \ 2k B T y ' J \2irk B T J y \ 2k B T K ' J V ' 

This local equilibrium distribution provides the correct averages for the first moments of /(r, v,w,i), this is, 

n(r,t)= / dvduf (r,v,LJ,t) 



n(r, i)u(r, t) — J dvdwv/o (r, v, w, t) 

n(r,t)f2(r,t) = J dvduiujf (r,v,u},t) (69) 

This, in turn, implies that 

J dwdujf 1 (r,w,u;,t)=0(V 2 ) 
J dvdwv/i(r,v,w,f) = C(V 2 ) 

J dvduujh (r, v, u), t) = C(V 2 ) (70) 

The procedure is now a bit different from the Chapman-Enskog method in Ref. because the inclusion of the 
spin variables produces new terms with different orders in gradients. We write the equation (|67f ) as follows 

dl 
4f 



^) fi-d t f 1 =d t f +v-Vfo-I[fo] (71) 

fo 



where we have neglected terms that are quadratic in gradients. We will check in the following that both sides of this 
equation are explicitly of first order in gradients. This linear equation ([7lJ) will be solved for /i and therefore we will 
obtain the solution of the FPBE (p6|) as fo + fi, up to terms of order 

We now consider each term of ( |71| ) separately. The temporal and spatial derivatives of fo can be computed to first 
order in gradients with the use of the balance equations ( p2[ ) and (|47]). Only terms of order V are to be retained, 
which amounts to use the balance equations with the averages of the quantities appearing in them evaluated with the 
local equilibrium ensemble. Therefore, we need to compute the local equilibrium average of the stress tensor II in 
the momentum balance equation, and the local equilibrium average of the two contribution in (p3^ ) to the equation 
for the angular velocity field. After using the molecular chaos assumption one easily obtains the following results 

II* - nk B TS^ 

n™ = - 7 «™ 2 ^ [At [cU^ + e^ a W] + B 2 [V-uV + + 8^]} + G(V 2 ) (72) 
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where we have defined 



A 2 = 



B 2 = 



D 



dRR 2 A(R) 



1 



D(D + 2) 



dKR 2 B(R) 



(73) 



The first contribution II^ produces an isotropic pressure term. Consistently with our assumption that the conser- 
vative forces are negligible this pressure is given by the ideal gas expression. The second contribution 11^° contains 
terms of first order in gradients. We arrange a bit this contribution by introducing the velocity gradient tensor 
(Vu)^ = c^u^ and its traceless symmetric and antisymmetric parts 



Vu + Vu J 



[Vu - Vu T ] 



1 

D 



-Vul 



(74) 



We have 



n 



no 



-777171 — 



A 2 {X7u A + + (A 2 + 2B 2 )Vu 5 + [A 2 + {D + 2)B 2 ] ^V-ul 



(75) 



The antisymmetric part of the total stress tensor in the local equilibrium approximation to first order in gradients is 
given by (as an axial vector) 



n 



AO 



2 A 2 



1 



jmn — — V x u f2 



(76) 



In a similar way one computes the quantities in (M) that appear in the balance equation for the spin (|47|). In 
particular, the last term in (|59| ) is also given by — 2lf T ^in the local equilibrium approximation to first order gradients. 

Substitution of the local equilibrium expressions for the stress tensor into the balance equations produce the Eulcr 
equations, 



dtn= -V?iu 

<9 f u = -(u-V)u - — -Vn + -V x 1 —n 2 Vt 

m n n 2 



d t fl = -(u-V)fi 4- 7yA 2 ri 



-X7 x u - n 

2 



(77) 



We have neglected a term of first order in gradients which produces a term of order V 2 in the momentum balance 
equation. We note that the time derivative of the angular velocity contains a term which is of zeroth order in gradients 
(the Q term in the last equation). 

With the help of the Eulcr equations and the chain rule, we can now compute the time and space derivatives of the 
local equilibrium distribution, to first order in gradient. The result is 



9*/o + v-V/ = /o 



Tfl I 

(v - u)(v - u) : Vu + (tu - 0)(v - u) : Vfl 



knT 



knT 



7n^ 2 ^-y — (v-u)-V x n^n + jAzn— (w-fl) 



-Vxu-O 
2 



(78) 



where the double dot : denotes double contraction. 

The next step is the calculation of I[fo\. To first order in gradients it is given by 



J[/o] = -fnA 2 



k H T 



fo 



-!-(v - u)-V x n 2 fl + (uj - SI)- 

2n 



ivxu- 



Therefore, after some happy cancelations the right hand side of (^Tj) has the simple form 



9t/o + v-V/ -/[/o] = /o 



Tfl I 

-V-u+ ^^r:(v - u)(v - u) : Vu + (co - fi)(v - u) : Vrj 



k B T 



ksT 



(79) 



(80) 
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which contains only terms of first order in gradients. 

Next, we consider the term dtfi in (fn]). We note that it is of first order in gradients due to the term of zeroth 
order in the Euler equation for the angular velocity, this is 



(81) 



where we have assumed that the dependence of /i on Q. appears in the combination u> — f2. This assumption will be 
confirmed a posteriori. 

The linearization of the functional /[/] might be easier to perform by expanding /[/0 + /1] to first order in gradients 
(at some point one uses Eqns. (|70|)). The final result for the left hand side of (|7l|) is 



I[fo + h] - I[M - d t fi = in 
where the operators are given by 

L 1 = — • v u 
ov 

OU! 



[A + B ]C T + ^C R 



fi + 0(V 2 



(82) 



k B T d 
m i9v 
k B T d 



(83) 



and the constants Bq are given by 



A = J dRA(R) 
Bo^^J dRB(R) 



Eqn. ([nj) can be written in compact form as 

Cfi = fa 
where the operator has the form 

C = 771 



TTi T 

V-u+ VV : Vu+ OV : Vf2 

kbT ksT 



[A + B ]C T + ^C R 



(84) 



(85) 



(86) 



and the peculiar velocities are V = v u, O = uj — fl. 

Equation (pq) is an inhomogeneous second order partial differential equation. In order to obtain a special solution 
of the inhomogeneous equations ( |85| ) , we introduce the following tensors 

1 



Dk B T 



T — o V 



With these quantities we write (pq) in the form 



Cfi = fo 



J : Vu +JS7-ul + T : Vfl 



(87) 



(88) 



The quantities tpTi) satisfy 
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^■ R fo"i^i> = 
£ T /o^ - -2/0.7 
>C fl /o^ = 

and therefore, a special solution of ( jS5[) is given by 



/1 = -/o 



1 



2 in {A + B ) 



J : Vu +JW-ul 



ln {A + B + f 



-T : Vf2 



(89) 



(90) 



as can be checked by substitution. 

Now it remains to obtain a general solution of the homogeneous equation Cfi = 0. The solution of this homogeneous 
equation is an arbitrary linear combination of fya, where a are the collisional invariants a = {l,v — u, uj — fi}. 
Nevertheless, the combination of and (p9f) imposes that the coefficients of the linear combination are zero. 



F. Transport coefficients 

The phenomenological theory of viscous flow of an isotropic fluid BTf] relates the trace tr [II] , the traceless symmetric 

II and antisymmetric U A parts of the stress tensor II with the linear velocity gradients and angular velocity in the 
following way 



tr[n] = -(V-u + p 



n 



-2r/Vu 



IT 4 = -2x]i 



-V x u - 
2 



(91) 



where the antisymmetric part is written as an axial vector. Here p is the isotropic hydrostatic pressure. The coefficients 
are the bulk viscosity £, the shear viscosity 7/ and the rotational viscosity r)R. 

We now compute the stress tensor (|5^) using the molecular chaos assumption (|6^) for the pair distribution function 
and the approximate solution fo + /1 for the single particle distribution function. This will produce II = IIo + III 
where the local equilibrium contribution IIo has been already computed in (|72|). Regarding the term III computed 
with /1 one observes that the only contribution which is of first order in gradients is Ylf which is computed along 
similar lines to Ref . |21f| . The final result is 



n 



K 



nk B T\ 



k B T =— s 
r Vu 



k B T 



<y[A + B ] 



Dj[A + Bo] 



Vul 



(92) 



The remaining contributions Ilf are of order V 2 and will be neglected. The final expression of the stress tensor in 
linear order of gradients is given by collecting fl79), (|76|) and (|9^) 



s tr[n] 



IT 



777m 
'jmn 



2D 



A 2 , (D + 2) 
2D 

^ + B 2 
2 2 



Bo 



1 D[A + B {) ] 



~/[A + B ] 



2 M 



1. 



-v x u - n 



Comparison of (^Tj) and ( |93| ) allows to identify the viscosities as 



Vu^ 



V-u + nk B T 



(93) 
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2 A 2 



>y[A + B }_ 



(94) 



In order to compare this expressions with those obtained by Marsh et al. in Ref. pT| , we should note that for the 
original DPD algorithm we have 



A(r) = 
B(r) = lu(R) 



(95) 



Simple substitution of (|95|) into (Q) shows that the transport coefficient (|94j ) coincide with those provided in Ref. 
ill. 



G. Transport equations 



Substitution of the stress tensor II = tr[II]l/Z) + n + IT A ( pi] ) into the momentum balance equation ( f42| ) produces 
the Navier-Stokes equations for a fluid with spin 27 (D = 3), 



p— u = -Vp + V-(2t7Vu s ) + V(C - 2?7/3)V-u + V x 
at 



2 m (n - ly x u) 



(96) 



where we have used the substantial derivative d/dt — <9 t +u-V. The last term in (|96j) is the gradient of the antisymmetric 
part of the stress tensor and describes the effect of the spin on the momentum transport. 

On the other hand by neglecting the term $ in the angular momentum balance equation ( |5l| ) [ p7| we obtain 



d t J = -V [Jv + r x n] 

which in combination with (^) produces the following balance equation for the spin density 

<9 t S = -V[Su] - 2U A 
which, implies the following dynamic equation for the angular velocity 



(97) 
(98) 



nl-fl = -2IT 4 

dt 



Substitution of Il A in (^lj) into this equation gives 



dt nl 



n — v x u 

2 



f2 V x u 

2 



(99) 



(100) 



The final closed set of equations for the hydrodynamic fields is given by Eqns. (|96|), ( J100] ), together with the continuity 
equation 



and the equation of state 



dt' 



p = ksTn 



c p 



(101) 



(102) 



where w e hav e introduced the speed of sound c. 

Eqn. (100) shows that the spin relaxes towards the vorticity with a relaxation time scale given by r = nl/Arjn p7| 
In the model of this paper, substitution of r]n in (EmI) gives the following time scale 



2"fnmA2 



(103) 
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H. Summary of kinetic theory 



In summary, it has been shown in this section that the macroscopic behavior of the fluid particle model i s hy dro d y - 
namical and the mass, momentum and angular momentum transport equations have been derived (Eqns. (101), (|9q), 
and ( |100| )). In this doing, explicit ex press ions for the transport coefficients in terms of the original model parameters 
have been obtained (Eqns. ( j94| ) and (103)). The equations here cited are the main results of the kinetic theory of the 
fluid particle model. 



VII. RESOLUTION ISSUES OF THE FLUID PARTICLE MODEL 



Within the picture of the Voronoi coarse-graining sketched in section II, it is possible to consider different levels of 
coarse- graining in which the number of atomic particles within a Voronoi cell is different. We expect that, provided 
that the number of atomic particles within the cell is large enough, the description of the hydrodynamic behavior will 
be more and more accurate as the number of Voronoi cells increases. In other words, we expect to reach a "continuum 
limit" as the number density of fluid particles goes to infinity. The discussion resembles that of the resolution 
in the numerical solution of partial differential equations. Actually, the resemblance can be made more accurate 
by comparing the structure of the equations of motion of the fluid particle model with those of smoothed particle 
dynamics. Smoothed particle dynamics is a Lagrangian discretization of the continuum equations of hydrodynamics 
that allows to interpret the nodes of the grid in terms of "smoothed particles" . 

For the case in which there is no coupling between the Navier-Stokes equation and the energy equation (the pressure 
does not depends on the temperature, for example), Takeda et al. H propose a discretization of the Navier-Stokes 
equations that produce equation of motions for the smoothed particles that corresponds exactly in structure with the 
postulated equations of motion of the fluid particle model in this paper. The correspondence is 



V(r) = 2 J\ 



"fA(r) 



W(r) 



1 






rjW 














mriQ 





W"(r) 



3 

W'(r) 



W'(r) 



(104) 



where, po, uq are the equilibrium pressure and number density, respectively, and W(r) is the weight function used in 
the discretization of the Navier-Stokes equation (the assumption that the density of all particles is almost constant 
has been taken). 

In this respect, the fluid particle model postulated in this paper is simply the smoothed particle dynamics with two 
additional bonus: 1) thermal noise is introduced consistently (that is, the fluid particle model can be interpreted as 
a Lagrangian discretization of the non-linear fluctuating hydrodynamic equations), and 2) the angular momentum is 
conserved exactly in the fluid particle model, in contrast with the smoothed particle dynamics model. The first bonus 
allows to apply smoothed particle dynamics to microhydrodynamic problems as those appearing in complex fluids 
where Brownian fluctuations are due to the fluctuating hydrodynamic environment. It can be also useful in studying 
the effect of thermal fluctuations near convective instabilities and, in general, in the study of non-equilibrium thermal 
fluctuations in hydrodynamic systems. The actual relevance of the second bonus will be discussed later. 

The comparison of SPD with the fluid particle model points out to an inconsistency that appears when usingsome 
particular selections for the weight function like the Lucy weight function § or a Gaussian weight function || . In 
these cases, it is easily seen that the function A(r) can become negative for certain values of r. This is unacceptable 
in view of Eqn. ([54j). From a physical point of view this means that if two particles are at a distance such that A(r) 
is negative, then the viscous forces will try to increase its relative velocities! 

In the derivation of the SPD model g], || it becomes apparent that the weight function W(r) must be normalized 
to unity, in order to have correct discrete (Monte-Carlo) approximations for integrals. If W{r) is normalized to 
unity, then one expects that by increasing the number density of smoothed particles one is increasing the numerical 
resolution of the simulation. The normalization implies that as the range of the weight function decreases with higher 
resolution, its height increases and, in the limit of infinite resolution it becomes a Dirac delta function W(r) — > S(r). 
Because the weight function is steeper when the resolution is higher, the time step used in the SPD model has to be 
reduced as the resolution increases. This is also encountered in any finite difference algorithm for solution of partial 
differential equations in order to maintain stability. 
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Let us investigate the effect of the resolution effects on the macroscopic parameters defining the fluid on hydrody- 
namic scales and which have been computed by means of the kinetic theory in the previous section. The parameters 
tha t ch aracterize the evolution of the velocity field are, as can be appreciated from (|9^), the speed of sound defined 
in ( 102 ) and the kinematic viscosities defined by v — n/p, Vb = C,/ p, an( l v r = Vr/ P- From ( pi]) they have the form 



u b = = 



7" 



2D 2D 2 



1 



jDn[A + B Q ] 



1 

v = — 
2 

vr = jn- 



C 7 nL4 + B ] 



(105) 



We are assuming, for the sake of the argument, that n — no, that is, the density field is constant. The conclusions, 
however, are valid in the compressible case also. 

Let us focus first on the dimensionless functions A(r), B{r) that determine the range of the dissipative and random 
forces. We expect that the clusters interact only with their neighbors, which are a typical distance A apart. Therefore, 
these functions will be of the form 
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(106) 



where a, b are functions that do not depend explicitly on A. This ensures that as the resolution is increased, the range 
of the force decreases, and this has the computationally appealing feature that the interaction between fluid particles 
remains always local. By using these scaling function and after a change to the dimensionless variable x = r/A, we 
have 
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where the dimensionless coefficients are given by 
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and do not depend on the resolution. By using (107) into (J105|) we obtain 
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We observe that all the dependence on the resolution (A or no) has been made explicit. In the limit of high resolution 
(A — * or no — > 00) the only contribution to the bulk if, an shear v viscosities comes from the kinetic contribution 
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that depends linearly on the temperature. This means that at zero temperature the system would not display any 
viscosity in the limit of high resolution. We find this behavior undesirable and we are lead to the conclusion that the 
friction coefficient 7 must depend on A. In particular, if we define 7 = 7A 2 (which has dimensions of a kinematic 
viscosity) and assume that 7 remains constant as the resolution varies, we will have 
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where the normalization ao = &o = 1 has been used as in the original DPD algorithm yjLA 
will take the form 

2m 702 



The relaxation time (102) 



(111) 



In this way, in the limit of high resolution (A — > 0) the viscosities are given essentially by 7 and the relaxation time 
goes to zero (note that the moment of inertia I must be of the form cx m x A 2 so t must tend to zero very fast). In 
the high resolution limit the spin becomes equal to the vorticity in a short time scale. The spin becomes a slaved 
variable and can be dropped from the description. Note also that in these situation the last term in the Navier-Stokes 
equation with spin ( |96| ) vanishes and one recovers the actual Navier-Stokes equation. This explains why, in SPD, the 
violation of angular momentum does not poses a serious problem for sufficiently high resolutions Q . If low resolutions 
are to be used in problems where the correct transfer of angular momentum is relevant (like in rotational diffusion of 
concentrated colloidal suspensions, for example), then the use of spin might suppose a real advantage. 

We have arrived at the conclusion that in order have a well-defined continuum limit the friction coefficient 7 must 
increase as the resolution increases. This can be understood physically in the following way. The number of particles 
in between two reference fluid particles at a given distance of each other increases as the resolution increases. If we 
require that the viscous interaction between these two reference particles must remain the same as the resolution 
increases, the mediating particles must interact stronger in order to transmit the same response between the two 
reference particles. From a mathematical point of view, the A 2 factor can be interpreted as the "lattice spacing" that 
is lacking in the original equations and that would be present in a numerical discretization of a second order derivative 
term. Preliminary simulation results for the DPD model (A(r) — 0) with energy conservation p9| shows that the 
correct continuum limit is obtained when the model parameter equivalent to 7 increases with A 2 jp6j. 

We would like to comment finally on an apparent inconsistency between SPD and the fluid particle model which 
is summarized as follows: if one discretizes the hydrodynamic equations on a set of points and then constructs the 
kinetic theory of these points, one would expect that the computed transport coefficients wou ld coincide with the input 
transport coefficients of the hydrodynamic equations. If one naively uses the results ( 104 ) in the calculation of the 
transport coefficients in ([m]), one arrives at an inconsistent result. The viscosities computed through the kinetic theory 
( |94| ) do not coincide with the input values. This could be traced back to the fact that the kinetic theory for the fluid 
particle model has been developed in the limit where no conservative forces are present, whereas the pressure ter m in 
the hydrodynamic equations (even for an ideal gas!) produces a conservative term given by the first equation in (104) 
in SPD. The kinetic theory with conservative forces is a bit more involved but the modifications can be summarized 
simply. The molecular chaos assumption ( |65| ) now involves the pair distribution function (which in the absence of 
conservative forces is equal to 1). This means that the parameters A2, B2 appearing in the transport coefficients will 
be modified by the presence of the pair distribution function within the integral defining these parameters. Also a 
new contribution to the transport coefficients arises due to the conservative forces. It is an open question whether 
these modified transport coefficients due to conservative forces do coincide with the input transport coefficients. The 
opposite case could also be possible simply due to the fact that the discretization procedure in SPD may induce 
"artificial viscosities" in the language of numerical resolution of the hydrodynamic partial differential equations. 

The fluid particle model is a consistent model by itself, without having to resort to the smoothed particle model 
for its validity. Actually, the fluid particle model, together with the kinetic theory developed in this paper has 
its advantages with respect to SPD: precise predictions can be made from the initial model parameters about the 
transport properties of the fluid. In this way, to obtain a prescribed fluid of known transport properties, one simply 
adjust the model parameters according to the formulae of kinetic theory (slight errors stemming from the failure of 
the molecular chaos assumption might play a minor role pl[). In SPD, on the contrary, the only way to specify the 
fluid is through the input transport coefficients in the original hydrodynamic equations. The discretization procedure 
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then produces a "fluid" whose transport properties do not in general correspond with those of the fluid intended to 
be modeled, and there is no systematic control on the appearance of artificial viscosities. 

It is apparent that this whole discussion can be applied to the DPD model, which is a particular case of the fluid 
particle model, and for which a kinetic theory has been formulated previously in Ref. |^T| . One of the main motivations 
for introducing shear forces between dissipative particles into the original algorithm of DPD was the identification 
of the following elementary motion between dissipative particles that produces no force in that algorithm. Let us 
focus on two neigbouring dissipative particles, the first one at rest at the origin and the second one orbiting in a 
circumference around the first one. This relative motion produces no force in DPD because the relative approaching 
velocity is exactly zero. Nevertheless, on simple physical grounds one expect that the motion of the second particle 
must drag in some way the first particle. This is taken into account through the shear forces in the fluid particle model 
presented in this paper. We note, however that this relative motion might produce a drag even in the original DPD 
algorithm if many DPD particles are involved simultaneously. The same is true for a purely conservative molecular 
dynamics simulation. The point is, of course, that the effect is already captured with a much smaller number of 
particles in the fluid particle model. 
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APPENDIX 



The derivation of the FPE is best achieved by considering the differential of an arbitrary function / to second order 



51 



At ^A 9 ^A d * ^A 9 f 

df = dvi h OVi - h duj 



d 2 f 



d 2 f 



3 a — 5 r dviduj-j- — - — 



A A d 2 / 
dwidvj- — - — 
ouJiOVj 



A A 9 V 
dujidujj- — - — 

OUJiOUJj 



(112) 



One then substitutes the SDE's ( |l9| ) and uses the Ito stochastic rules ( |L7| ) keeping terms up to order dt (the cross 
terms involving positions have been neglected in (112) on account of the fact that dr is already of order dt). Then 
after averaging with respect to the distribution function p(r, v,w,t), one performs a partial integration and uses the 
fact that / is arbitrary, to obtain the Fokker-Planck equation in the form 



d t p{r, v, oj; t) = [L c + L T + L R ] p(r, v, oj; t) 
where we have defined the operators 
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(114) 



The operator L c is the usual Liouville operator of a Hamiltonian system interacting with conservative forces F c . We 
need to arrange a bit the operators L T and L R by using the Ito rules (18) 
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The second order tensor T^V- , satisfies 
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If we define 
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then the following identities are obtained 
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By using these results into (114) the operators take the following compact form 
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where we have introduced the vector operators 
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